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ABSTRACT 

NGC 604 is the second most massive H n region in the Local Group, thus an important laboratory for 
massive star formation. Using a combination of observational and analytical tools that include Spitzer 
spectroscopy, Herschel photometry, Chandra imaging, and Bayesian Spectral Energy Distribution 
fitting, we investigate the physical conditions in NGC 604, and quantify the amount of massive star 
formation currently taking place. We derive an average age of 4 ± 1 Myr and a total stellar mass of 
1-61i'q x 10 5 Mq for the entire region, in agreement with previous optical studies. Across the region we 
find an effect of the X-ray field on both the abundance of aromatic molecules and the [Si n] emission. 
Within NGC 604 we identify several individual bright infrared sources with diameters of about 15 pc 
and luminosity weighted masses between 10 3 Mq and 10 4 Mq. Their spectral properties indicate that 
some of these sources are embedded clusters in process of formation, which together account for ^8% 
of the total stellar mass in the NGC 604 system. The variations of the radiation field strength across 
NGC 604 are consistent with a sequential star formation scenario, with at least two bursts in the last 
few million years. Our results indicate that massive star formation in NGC 604 is still ongoing, likely 
triggered by the earlier bursts. 

Subject headings: H n regions - infrared: ISM - ISM: individual: NGC 604 - stars: formation 



1. INTRODUCTION 

Our current understanding of massive star-forming re- 
gions remains poor, despite their importance in the struc- 
ture and evolution of galactic systems, due to their strong 
feedback. From an observational point of view, a num- 
ber of reasons make the study of massive star formation 
a challenging topic in astrophysics, a s pointed out in the 
compr ehensive review on the issue by lZinnecker &: Yorkd 
(2007). Regions of massive star formation are highly em- 
bedded in thick layers of dust during the crucial early 
stages of their existence. In addition, these early stages 
are short-lived, leaving little time for the study of their 
evolution. Another problem is the lack of spatial reso- 
lution in the observations. Most of extragalactic giant 
star forming regions are often contained within a single 
pixel. Finally, most massive stars form in close proxim- 
ity to each other, and their mutual influence via gravi- 
tational interaction, powerful outflows, supernova events 
and strong winds contributes to the complexity of the 
problem. 

Observations of nearby star forming regions provide 
an excellent laboratory for the study of a particular as- 
pect of this interaction, namely the triggering of new 
star formation events by gas compression resulting ei- 
ther from a ionization shock front created by the radi- 
ation field of a previous generation of stars, or by the 
supersoni c shocks of a supernova ev ent. In the classical 
theory by Elmcgrc en &: Ladal (|1977l ). an ionization front 
compresses an adjacent layer of molecular gas and heats 
it, producing a gravitational instability that will even- 
tually result in a new generation of stars. In order to 
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test this and other theories of triggered star formation, 
it is important to unequivocally determine and quantify 
the amount of currently ongoing star formation in Giant 
Hn Regions (GHnRs), and its relation to the ionizing 
radiation from previous generations of stars. 

After 30 Doradus, NGC 604 is the second most mas- 
sive GH nR in the Local Group. Locate d in the Triangle 
Gala xy (M33) at a distance of 0.84 Mpc (jFreedman et al.l 
1991), it harbors several associations of massive stars dis- 
tributed across an area of about 200 pc on the side, 
dominated by a c luster containing ~ 200 OB stars 
(jHunter et al.lll996l . and references therein). These as- 
sociations have excavated a complex system of filaments 
and cavities of ionized material surrounded by photon- 
dominat ed regions (PDRs) and mo lecular gas (see, for 
example Rela no &: KennicuttJ 120091 for a discussion on 
the spatial distribution of emission in the region). Opti- 
cal stud ies reveal an age o f the region between 3 and 
5 My r (jHunter et al.l 119961 : iGonzalez Delgado fc Pered 
2000) and a total s t ellar mass of (3.8 ± 0.6) x 10 5 M 
(jEldridge fc Relanol 120111 ). Individual CO molecular 
clouds have been detected with sizes between 5 and 
29 pc and with masses of between 0.8 x 10 5 Mq and 
7.4 x 10 5 Mq (jMiura et al J 120101) . Measured values of 
the average optical extinction in the r egion vary from 
Ay = 0.3 (iRelano fc Kennicuttl 120091) to A v = 0.5 
(jChurchwell fc Gosslll999r i. iRelano fc Kennicuttl (2009) 
also derive a star formation of about (5.7±0.4) x 10 5 Mq , 
which can be interpreted in terms of the total mass of a 
coeval stellar population. 

Several attempts have been made to quantify the 
total amount of ongoing massive star formation in 
NGC 60 4. Using near-infr ared (NIR) color-color di- 
agrams, iBarba et al.l (|2009l ) identify several Massive 
Young Stellar Object (MYSO) candidates, with well- 
defined infrared excesses. They also note that these 
candidates coincide spatially with radio-peak structures, 
reinforcing their star-forming nature. More recently, 
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I Farina et~a l. (2012) report the discovery of sources with 
near infrared excess within the infrared-bright ridges 
surround ing the ionized gas, and a ssociate them with 
MYSOs. iRelano fc Kennicuttl (|2009T l argue that the red- 
dening observed towards some prominent sources in the 
region can be explained by the existence of foreground 
molecular material, but they do not rule out the possi- 
bility of those sources being embedded sites of star for- 
mation. The use of infrared spectrometry combined with 
physical modeling provides a powerful tool to distinguish 
between foreground extinction and embedded star for- 
mation, by fitting the spectra with independent models 
that correspond to each situation, or to a combination of 
both. This complements photometric methods that use 
near-infrared color-color and color-magnitude diagrams 
to discri minate between the tw o cases, as the one de- 
scribed in lLada fc Adams! (|1992ft . 

In this paper we perform a comprehensive analysis of 
the physical conditions in NGC 604 in the context of its 
evolutionary status, and investigate the presence of ongo- 
ing massive star formation in the region. We use infrared 
spectral and photometric data from the Spitzer and Her- 
schel Space Telescopes, complemented with archival X- 
ray and optical data, as well as a set of analytical tools 
to interpret them. We report the discovery of individ- 
ual infrared knots whose derived masses are consistent 
with them being stellar cluster in process of formation. 
This is also supported by Bayesian fitting of the Spectral 
Energy Distribution (SED) of the region, which points 
to the presence of a significant component of embedded 
objects in the region. We derive line emission and contin- 
uum maps of NGC 604 and use them to assess the varia- 
tions in radiation field strength, ionization levels and ex- 
tinction across the GH nR, and find that our results are 
consistent with a sequential star formation history in the 
last ~ 4 Myr. We also discuss some additional findings 
regarding the role of X-ray emission in the enhancement 
of both [Sin] atomic emission and 17 /im PAH emission. 

The paper is structured as follows. In $2]we present the 
IRS data and describe the data reduction to obtain the 
maps and the spatially integrated spectrum of NGC 604. 
Sj3] describes the resulting maps and spectra, as well as 
the tools used to extract physical information from the 
observations. In S|4]we discuss the results of our analysis 
in terms of the current evolutionary stage of NGC 604 
and the presence of ongoing star formation in the region. 
We conclude with a summary of our main results in Sj5] 

2. DATA REDUCTION AND ANCILLARY DATASETS 

Most of our analysi s will be based on Spitzer Infrared 
Spectrograph (IRS) (|Houck et all I2004D data of a re- 
gion encompassing the bulk of infrared emission from 
NGC 604. However, to provide more robust constraints 
on the physics of the region, we also use complemen- 
tary pho tometry from the Spitzer Infrared Array Camera 
(IRAC) (|Fazio et al.ll2004D and from the Herschel Space 
Observatory Phot odetector Array Cam era and Spec- 
trometer (PACS) (|Poghtsch et all 12010ft . Additionally, 
we use archival images from Hubble Space Telescope Wide 
Field and Planetary Camera 2 (WFPC2) and Chandra 
X-ray Observatory Advanced CCD Imaging Spectrometer 
(ACIS). 

2.1. IRS data 



The Spitzer IRS provides unprecedented spatial resolu- 
tion and sensitivity as compared to any previous mid-IR 
spectroscopic observations of the NGC 604 region. The 
spectral resolving power ranges from R ~ 60 at short 
wavelengths to R ~ 120 at the long- wavelength edge. 
With a 512 s exposure, point source sensitivity limit be- 
tween 0.1 and 10 mJy in the 5-38 fim range, the IRS is 
about 100 times more sensitive than the Infrared Space 
Observatory (ISO), while the spatial resolution is a factor 
of 10 larger. 

All the observations presented in this section were ob- 
tained in January 2006 using the IRS mapping and star- 
ing modes for spectroscopy, as part of the program Com- 
parative Study of Galactic and Extragalactic Hu Regions 
(P. I. J. Houck). The mapping mode consists of the ac- 
quisition of slit spectra using a grid of positions around 
a central target. Only the low resolution modules short- 
low (SL1, SL2) and long-low (LL1, LL2) of the IRS were 
used for the NGC 604 map presented in this paper. For 
the SL modules, 12 slit pointings were made with each of 
the two spectrometer orders covering an area on the sky 
of about 55" x 40", which corresponds to a physical scale 
of about 225 pc x 160 pc at the distance of NGC 604 (see 
Fig- HI)- The slice width for the SL modules is 3". 6, cor- 
responding to a pixel scale of 1".85 px _1 (7.5 pc px _1 ). 
For the LL modules, the slice width is 10". 5, which cor- 
responds to a pixel scale of 5".08px _1 (20.5pcpx _1 ) and 
6 pointings were made with each order to cover an area 
of about 720 pc x 205 pc. 

In addition to the spectral map with the IRS low reso- 
lution modules (lores data hereafter) , staring mode spec- 
tra were also obtained with the IRS for three specific lo- 
cations within the region (hires data here after). In this 
"point and shoot" mode, the targets are placed in the 
center of one or several of the slits for a specified inte- 
gration time. In staring mode, only the high resolution 
modules of the spectrometer were used, namely short- 
high (SH) and long- high (LH). The wavelength coverage 
is shorter (between 9.9 /im and 38.0 /im), but the re- 
solving power is significantly higher (i? 600). The lo- 
cations of the staring mode observations are within the 
area of the spectral map and correspond approximately 
to the positions of the peaks of 8 fim emission, associated 
with PAH emission, as discussed later. The high resolu- 
tion slits are wider than their low resolution counterparts 
(4". 7 and 11". 1 for SH and LH, respectively); therefore 
their spectral apertures are also larger than the pixel 
scale of the lore spectral map. 

The area of the lores map and the locations of the 
staring mode hires spectra within the map are shown in 
the right panel of Fig. [2] 

2.1.1. Extraction of the spectra 

For the staring mode data, the Spitzer Science Center 
Tool IRSCLEAN was used to remove cosmic rays. Then, 
SMART v. 8.0 was used with the full aperture mode for 
extended sources, based on the size of the slits compared 
to the source NGC 604, which is large enough to be con- 
sidered extended. The observation cycles were co-added 
and the sky background removed using an additional off- 
source background exposure taken as part of the cam- 
paign. We have applied a scaling factor to the fluxes 
extracted in the LH module to match the overlap region 
in the SH module. This is an aperture correction to ac- 
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count for the larger size of the LH slit with respect to the 
SH slit. The scaling factors for the LH spectra towards 
fields A, B, and C are 0.28, 0.41, and 0.41, respectively 
(see Fig. [2]). 

For the spectral map, we extracted the spatially inte- 
grated spectrum of the region for all four IRS orders using 
a extraction aperture corresponding to the SL coverage 
of the map (see Fig. [2J , which includes the bulk of the IR 
emission of NGC 604. Some of the extended filamentary 
structure of the region, which is one order of magnitude 
dimmer than the peaks of emission, is left outside this 
area. The low surface brightness of these filaments and 
the bubble-like geometry that we will assume when mod- 
cling the region implies that this will be unimportant for 
the purposes of the present analysis. The background 
subtraction was performed using the order of the SL and 
LL modules that was not centered at the source during 
the corresponding exposure. CUBISM performs the data 
cube build-up once the background and the correct slit 
pointings have been provided. The background levels are 
not exactly the same for the SL and LL modules. The 
reason is that the orientation of opposite orders of the slit 
is different for different modules. In the case of the low 
resolution modules, the background picks up some emis- 
sion from M33's spiral arm (See Fig.[T]). This background 
emission from the spiral arm is our best estimate for the 
sky levels in the SL modules. For reference, the measured 
IRAC 8 /im flux levels on the spiral arm (~ lOMJy sr _1 ) 
are only 20% higher than the sky level outside the arm 
(~ 8MJy sr _1 ), and correspond to about 7% of the peak 
of 8 /im emission in the region. 



TABLE 1 
Integrated IRAC photometry 
of NGC 604 




Fig. 1. — The target and background SL slit pointings super- 
imposed on the IRAC 8 fira image of NGC 604. Color coded 
are the IRAC fluxes from about ~ 8 MJy sr — 1 (dark blue) to 
~ 120 MJy sr _1 (yellow). Emission from the spiral arm of M33 
is visible towards the south-east of the Hll region, in the region 
where the background was taken. See text for a discussion on the 
flux levels. 

Spectra can be extracted from each spatial pixel of the 
resulting data cube. The spatially integrated spectrum 
over the entire aperture is obtained by summing up the 
individual spectra of all resolution elements. We use the 
same extraction area for all modules and orders using the 
tools provided by CUBISM to make sure that the final 
spectrum for each order corresponds to the same physical 
region. 

In addition to the spectrum of the integrated region, 
we have extracted the lores spectra of individual sources 



Wavelength [/tm] 


Flux [Jy] 


3.6 


0.067(0.002) 


4.5 


0.062(0.002) 


5.8 


0.322(0.010) 


8.0 


0.922(0.028) 



that we discuss in §3.1.21 For this extraction we use an 
aperture of a single SL pixel (T'.85). The CUBISM soft- 
ware allows us to extract the spectra in the LL spectral 
module with this small aperture, by scaling down the 
fluxes of the larger LL pixels to the SL sizes. This is 
equivalent to applying a scaling factor to match the dif- 
ferent orders of IRS spectra, and introduces additional 
flux uncertainties in the long wavelength modules, where 
the PSF is not well sampled. Nonetheless, we consider 
this scaling a good approximation of the actual fluxes 
at smaller scales, since the SL pixels are not sufficiently 
small to resolve the sources, and the LL pixels are not 
large enough to include more than one source. 

By selecting a single pixel aperture, we have cho- 
sen to loose some spatial information on the individ- 
ual sources, because a single pixel samples only half of 
the PSF FWHM, and in exchange we avoid off-source 
flux contamination. More important than a full sam- 
pling of the PSF for our analysis, are the variations of 
the PSF with wavelength, that might introduce arti- 
facts in th e measured spectral fea t ures. Using calibra- 
tion data, iPereira- Santaell a et al.l ((2010) characterized 
the PSF variations with wavelength for the IRS recon- 
structed PSFs. Apart from an undulating behavior of the 
PSF size with wavelength, that they attribute to align- 
ment issues and to the reconstruction algorithm used, 
they compute variations of less that 10% on the PSF 
FWHM for the SL module. Such variations are below 
the observational errors in our data, described in £ 12 .31 
and I3T2"1 The PSF of the LL module is affected by 
fringing and is difficult to characterize. 

2.2. IRAC photometry 

As part of the same Spitzer program, IRAC maps of 
the NGC 604 region were obtained at 3.6 /im, 4.5 /zm, 
5.8 /zm, and 8 /im. The pixel scale for these maps is 
1.2" px _1 , and they cover an area of about 5' x 9', sev- 
eral times larger than the area covered by the spectral 
map. For the purpose of this paper we have extracted the 
integrated flux of each map within a rectangular aperture 
area equal to the extraction area of the IRS spectral map. 
We perform this extraction using the FUNTOOLS pack- 
age for the SAO ds9 software. The sky background is 
estimated from the map by measuring the flux in a box 
of the same size as the map, but shifted to the west, to an 
area where no source emission is observed. In Table [T] we 
list the measured fluxes. The listed uncertainties corre- 
spond to absolute flux calibration uncertainties (~ 3%), 
which are derived for point sources taking several sys- 
tematic effects into account, such as array position de- 
pendence, pixel phase dependence, color correction, an d 
aperture correction, as described in lReach et al.l {2005). 

The IRAC maps provide a sharper view of the region 
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as compared to the IRS spectral map at specific wave- 
lengths, and will allow the identification of interesting 
sources. 

2.3. PACS photometry 

Imaging maps of the host galaxy M33 have been ob- 
tained with the Herschel PACS instrument using the 
green (100 /zm) and red (160 /im) filters as part of the 
HERM33ES Herschel key project, and comprehensive de- 
scription of_thejdjia_redu^tion and obained maps can be 
found in lBoauien et al.l (|2011[ ). The maps were obtained 
with a slow scan speed of 20" s -1 , and cover a total area 
of about 70' x 70' . Here we use the integrated photometry 
of an area of the PACS maps equivalent to the size of the 
IRS spectral maps (right panel of Fig. [2]). The pixel sizes 
are 3". 2 for the green band and 6". 4 for the red band; 
about 2 and 4 times the IRS-SL pixel size. The obtained 
fluxes, integrated over the entire area of the map, are 
JWn =39.7±4.0JyandF 160Aim = 30.1 ±3.0 Jy. The 
10% uncertainty comes from a combination of absolute 
calibration errors, uncertainties associated with differ- 
ences in the PSF at 100 /an and 160 /im, and different 
pixel sizes in the two bands that lead to aperture un- 
certainties. The rms noise levels of the PACS maps are 
2.6 mjy px~ 2 and 6.9 mJy px -2 . 

2.4. HST-WFPC2 F555W data 

We use the optical images obtained at 0.55 /im with 
the Hubble WFPC2 using the F555W filter, described in 
iHunter et al.l (|1996|) . At angular resolutions of 0".l, this 
optical map reveals the location of the massive ionizing 
clusters that provide the radiative input for the NGC 604 
system. 

2.5. Chandra X-ray Observatory- ACIS data 

We use archival ACIS X-ray data. The data were 
taken as part of the Chandra proposal "The Giant Ex- 
tragalactic Star-Forming Region NGC 604" (Proposal ID 
02600453, P.I. F.Damiani), and consist of a soft (0.5- 
1.2 keV) X-ray image of the entire nebula, with an expo- 
sures time of 90 ks. The pixel scale is 1" px _1 . 

3. ANALYSIS 

We analyze the multi-wavelength observations de- 
scribed above using a set of analytical and statistical 
tools that are based on physical models of the region, 
and that we will describe shortly. The models compute 
the radiative transfer of the UV radiation as it traverses 
the ionized gas and molecular material around the H n 
region. They also compute the dynamical evolution of its 
the expanding Hn region. We use these tools to derive 
physical properties of the region, such as dust tempera- 
tures, total stellar mass, hardness of the radiation field, 
and ionization state of the gas. In this section we present 
the obtained maps and spectra and describe the analyt- 
ical tools that we use. 

3.1. Distribution of the emission 
3.1.1. Overall distribution 



0.55 /zm (blue) . The blue channel shows the photo- 
spheric emission from young massive stars, most of which 
belong to the dense clu s tered structure labeled as "clus- 
ter A" by Hu nter et al.l (|1996l ). Other, more spread stel- 
lar associations are seen next to the bright lobes of in- 
frared emission, most notably within fields B and C. 
The 8 /zm emission traces warm dust and the 7.7 /im 
PAH feature, and hence the location of the PDRs, while 
the 3.6 /zm traces the 3.3 /zm PAH feature and pho- 
tospheric emission fro m stellar populations older than 
lOMyr, mostly giants (jHelou et all 12004 ICannon et all 
2006). The filamentary and shell- like structure of the 
region observed at optical wavelengths is reproduced in 
t he mid-infrared. 

iRelaho fe Kennicu tt (2009) have shown that the 8 /zm 
emission delineates the Ha shells in the boundaries be- 
tween cavities. Most of the infrared emission in our 
maps comes from two lobes oriented in a SE-NW di- 
rection and that constitute the bulk of the emission at 
the IRAC bands. These lobes coincide with the edges of 
two main cavities observed i n the region (cavities Bl and 
B2 in lTiillmann et all [2008). with diameters of approxi- 
mately 50 pc each. They also coinci de with the position 
of bri ght radio knots identified by IChurchwell &; Gosi 
(1999). Within these cavities sits the majority of the 
luminous stars observed by HST (blue stellar sources in 
Fig. A n exhaustive X-ray stu dy of the region car- 
ried out bv lTullmann et all (|2008f ) showed that the dis- 
tribution of high energy photons inside these cavities is 
consistent with them being shaped by the mass loss and 
radiative pressure of about 200 OB stars. The eastern 
portion of NGC 604, on the other hand, seems to be an 
older part of the system, with X-ray emission consistent 
with a more e yolved population. Based o n optical spec- 
troscopy, also iTenorio-Tagle et al.l ()2000D report a kine- 
matical difference between the eastern and western parts 
of NGC 604 and claim that the eastern part is not af- 
fected by the ionized cluster. 

The two lobes have similar surface brightness at 8 fj.ui 
but they show differences in their morphology. At the 
IRAC wavelengths, the SE lobe shows a series of sub- 
condensations, some of which are notably brighter at 
3.6 /im while the NW lobe shows a more uniform distri- 
bution of emission. CO maps of the region reveal several 
molecular clouds in the region, with a CO-bright cloud 
peaking near the SE infrared lobe and extending south- 
wards, an d a dimmer and smaller cloud peaking near the 
NW lobe (I Wilson fc Scovillelll99l . A number of MYSO 
candidates have been identified along these two lobes as 
sourc es with NIR excess (see, for example iFarina et al.l 
12012ft . By combining observational and SED modelling 
tools, in this chapter we will provide additional evidence 
that supports the star formation scenario, and quantifies 
its contribution to the total stellar mass in the region. 

3.1.2. Individual Sources. 



In Fig. 2(a) we show a three-color map of NGC 604 
composed from the 8.- /im (red) and 3.6 /im (green) 
IRAC channels together with the WFPC2 image at 



Fig. 2(b) shows a three color map of NGC 604 with the 
main extraction areas indicated (see Table [5]). The ex- 
traction area for the spectral map includes the bulk of the 
emission in the IRAC bands and the stellar cluster. We 
have selected 7 individual sources of interest in the spec- 
tral map, including the well defined sub-condensations 
that are visible in the IRAC images. These sources are 
indicated in Fig. 2(b) with their positions listed in Table 
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TABLE 2 
IRS targets in NGC 604 
Lores sources and hires fields 



Source 


RA 


Dec 


Src 1 


Ih34m33.43s 


+30°46'48.50" 


Src 2 


Ih34m33.67s 


+30°46'51.01" 


Src 3 


Ih34m33.70s 


+30°46'54.86" 


Src 4 


Ih34m33.43s 


+30°46'57.11" 


Src 5 


Ih34m33.56s 


+30°47'03.00" 


Src 6 


Ih34m32.73s 


+30°47'09.32" 


Src 7 


Ih34m31.99s 


+30°46'59.95" 


A 


Ih34m33.6s 


+30°46'51" 


B 


Ih34m33.6s 


+30° 46' 58" 


C 


Ih34m32.4s 


+30°46'59" 



Fig. 2. — (a) Three-color Hubble/Spitzer map of NGC 604. The IRAC 8 (im band traces the PDR material, while the 3.6 /im band 
traces both PAH and photospheric emission from young stars. The WFPC2 image shows the location of the hot massive stars, (b) the 
same map with superimposed apertures of the IRS slits (A, B, and C) and labels for seven sources of interest based on their IRAC colors. 
We discuss these sources in the text and indicate their location in the map. The yellow rectangle corresponds to the extraction window for 
the IRS map and for the photometry, while the smaller magenta boxes show the field location of the IRS SH slit for the high resolution 
observations. The LH slits had the same locations as the SH slits. North is up, east is to the left. 

3.2. Infrared Spectra 
3.2.1. Extracted Spectra 

In Fig. [3] we show the lores IRS spectra of the inte- 
grated region and the selected sources of Table [5J The 
vertical axis is in units of flux density (uFu), and the 
spectra have been scaled by progressive factors of 10 to 
show the direct comparison. The integrated spectrum of 
NGC 604 has prominent PAH emission and some of the 
nebular lines detected are from species such as [Ann], 
[Ne III] , [Ne n] , [S iv] , [S in] , and [Si n] . The thermal con- 
tinuum increases monotonically in the IRS range. 

Fig. @] shows the high resolution spectra of the star- 
ing mode targets. The wavelength coverage of the hires 
modules is smaller than the lores case, from ~ 10 /im to 
~ 37 /im, but the higher spectral resolving power allows 
the resolution of lines not seen in the lores spectra, such 
as the lines resulting from pure rotational transitions of 
molecular hydrogen, H 2 S(2) at 12.3 [im and H 2 S(l) at 
17.0 /zm. These lines are usually hard to detect on top 
of a strong continuum, due to the fact that they arise 
from quadrupolar rotational transitions, which are in- 
trinsically weak. 

3.2.2. Continuum emission 

In order to characterize the spectral slope of the ther- 
mal continuum, we measure the flux densities at 15 /jm 
and 30 /jm using a range of wavelengths around the cor- 
responding central wavelength containing about 20 res- 
olution elements (14.75 /im -15.25 /im for the 15 /im 
measurement and 29.5 /im -30.5 /zm for the 30 /im mea- 
surement). We calculate the ratio Fi 5fim /F 30fim for each 
source as well as the integrated spectrum. The spectral 
slopes give an indication of the dust temperature. Higher 
values of Fupm/Faoftm are associated with a hotter com- 
ponent of the dust, whereas lower values of this ratio 
indicate colder dust temperatures. We list the measured 
fluxes and slopes in Table [31 We have used a nominal 
10% uncertainty in the IRS fluxes, which accounts for 



[2 Fields A and B of the hires observations include some 
of the selected sources of the lores map. Specifically, field 
A includes sources 1 and 2, while field B includes source 
4. Field C does not include any of the selected lores tar- 
gets, but belongs to the NW infrared lobe, and is close to 
source 7. Besides, fields A, B, and C correspond to radio 
fields B, A, and C, respectively, of IChurchwell fc Gossl 
Jl99flh. while source 6 is slightly shifted from the cluster 
A of iHunter etHl (fl996h . 

We have fitted Gaussian profiles to the flux distribution 
of the observed sub-condensations to investigate their 
spatial extension due to the PAH emitting regions. Our 
fits reveal that they have projected diameters of about 
3". 6 at 8 /zm, as measured from the FWHM of the Gaus- 
sian fits. The diffraction-limited resolution of the IRAC 
camera at this wavelength is 1".71, and hence we con- 
clude that the sub-condensations are spatially resolved at 
8 /im, having a diameter of at least two IRAC resolution 
elements. At the distance of NGC 604, their projected 
sizes correspond to physical diameters of about 15 pc; 
hence, comparab le to the size of a typical galactic giant 
molecular cloud (jFukui fc Kawamurall2010D . 
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Fig. 3. — IRS spectra of the integrated NGC 604 region and 
the individual sources listed in Table [2] The position of some 
prominent mid-infrared fine structure lines are indicated, as well as 
the location of the PAH bands. The small feature seen at ~20 fim 
in sources 3 and 5 is an artifact from the data reduction. 



TABLE 3 

SED CONTINUUM SLOPE FOR ALL SEVEN SOURCES AND 
THE INTEGRATED MAP. 



Source 




F30Aim 






[MJy sr" 1 ] 


[MJy sr" 1 ] 




Src 1 


102.8(10.3) 


432.1(43.2) 


0.238(0.034) 


Src 2 


65.4(6.5) 


310.4(31.0) 


0.211(0.030) 


Src 3 


95.8(9.6) 


582.8(58.3) 


0.164(0.023) 


Src 4 


95.9(9.6) 


447.5(44.8) 


0.214(0.030) 


Src 5 


25.8(2.6) 


258.5(25.9) 


0.100(0.014) 


Src 6 


10.1(1.0) 


99.2(9.9) 


0.101(0.014) 


Src 7 


113.5(11.4) 


493.9(49.4) 


0.230(0.033) 


NGC 604 


16.8(1.7) 


100.9(10.1) 


0.166(0.023) 



the absolute and relative calibration of the instrument, 
as well as systematic errors due to specific observing con- 
ditions. 

3.2.3. PAH emission 

We have measured the strengths of the individua l 
PAH features using the PAHFIT tool (|Smith et al.ll2007D . 
which decomposes the IRS spectrum in individual contri- 
butions from dust thermal continuum, PAH features, and 
fine-structure lines, as well as rotational lines of molecu- 
lar hydrogen. In TableHlwe list the measured continuum- 
subtracted strengths of the different PAH components, 
and their uncertainties as derived with PAHFIT. The 
7.7 /xm strength was obtained by adding the PAHFIT 
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Fig. 4. — High resolution spectra of the three fields labelled A, B, 
and C in Fig. [2] The prominent nebular and molecular hydrogen 
lines are indicated. 

7.4 /xm, 7.6 /im, and 7.8 /im features. Similar combina- 
tions were performed for the 11.3 /im feature (using the 
11.2 /im and 11.3 /im strengths) and for the 17 /im fea- 
ture (using the 16.4 /xm, 17.0 /xm, 17.4 /xm, and 17.9 /xm 
strengths). 

In general, the ratios between PAH feature strengths 
have similar values for all sources, with source-to-source 
variations of between 6% and 14%. An interesting 
exception is the ratio of the PAHi7 Mm feature to the 
sum of all the other PAH features, EPAH (EPAH = 
PAH 6 . 2M m + PAH 7 . 7M m + PAH 8 .6 M m + PAHu.3^). This 
ratio shows source-to-source variations of about 60% and 
peaks strongly towards source 3, which is also a local 
peak of soft X-ray emission. In Table [5] we list the 
PAHi 7Aim /SPAH ratio and the Chandra- ACIS soft X-ray 
fluxes for our sources. For the ACIS fluxes we assume a 
nominal 20% uncertaint y, based on Chandra catalogues 
of point sources such as IServillat" 

3.2.4. Nebular line emission 

Fine-structure emission lines are an important diag- 
nostic of the physical conditions in star-forming regions. 
Their strengths and ratios constrain physical parame- 
ters such as the hardness of the radiation field, gas 
density, or the presence of shocked gas. In the mid- 
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TABLE 4 

PAH FEATURES STRENGTHS 



Source 


PA LT,. ~ 


PA TT™ ~ 


PA W~ r- 








10- 6 x[W mr 2 sr- 1 ] 


10- 6 x[W m" 2 sr- 1 ] 


10" 6 x[W m~ 2 sr- 1 ] 


10- 6 x[W in- 2 sr" 1 ] 


10~ 6 x[W m- 2 sr- 1 ] 


Src 1 


2.56(0.10) 


7.96(0.54) 


1.28(0.12) 


1.33(0.16) 


0.85(0.26) 


Src 2 


1.76(0.07) 


5.49(0.38) 


0.87(0.09) 


1.07(0.11) 


0.73(0.17) 


Src 3 


1.24(0.07) 


2.82(0.33) 


0.51(0.09) 


0.83(0.14) 


1.15(0.24) 


Src 4 


1.51(0.07) 


4.25(0.35) 


0.64(0.09) 


0.76(0.13) 


0.88(0.23) 


Src 5 


0.93(0.03) 


2.42(0.19) 


0.50(0.04) 


0.52(0.05) 


0.24(0.07) 


Src 6 


0.14(0.01) 


0.48(0.03) 


0.06(0.01) 


0.11(0.01) 


0.06(0.03) 


Src 7 


1.64(0.08) 


5.14(0.41) 


0.91(0.10) 


1.04(0.14) 


0.26(0.07) 


NGC 604 


0.40(0.02) 


1.16(0.09) 


0.19(0.02) 


0.25(0.03) 


0.11(0.04) 



TABLE 5 

Relation between the ratio of PAH 17 /^m to all the 

OTHER PAH FEATURES AND SOFT X-RAY EMISSION. 



Source 


PAH 17flm /(£PAH) a 


Fx aoft [10" 9 erg cm" 2 s" 1 ] 13 


Src 1 


0.065(0.020) 


0.84(0.17) 


Src 2 


0.079(0.019) 


1.37(0.27) 


Src 3 


0.213(0.047) 


3.70(0.74) 


Src 4 


0.123(0.033) 


2.50(0.50) 


Src 5 


0.055(0.016) 


1.64(0.33) 


Src 6 


0.076(0.038) 


1.65(0.33) 


Src 7 


0.030(0.002) 


0.91(0.18) 


NGC 604 


0.055(0.020) 





a Fluxes are from PAHFIT. EPAH is obtained as the sum of the 
PAH features at 6.2 /im, 7.7 ^m, 8.5 (im, and 11.3 fJ.ni. 
b Measured from Chandra-ACIS data. 



IR, several forbidden lines are observable from species 
such as [Ne n] , [Ne in] , [S m] , and [S iv] . In Table EO we 
list the continuum-subtracted nebular line strengths in 
NGC 604, as measured with the PAHFIT tool, as well as 
their uncertainties from the fit, also returned by PAH- 
FIT. 

We have also measured the continuum-subtracted line 
strengths in the hires data for targets A, B, and C in the 
right panel of Fig. [5] We have fitted Gaussian profiles 
to the lines, this time using the built-in tool for that 
purpose included in the SMART software. We show the 
resulting line strengths in Table [7] 

Fig. [5] shows continuum-subtracted maps of the 
[Siv]10.5/im and [Ne ii]12.8/im fine-structure lines. 
Within the maps spatial resolution there is spatial coin- 
cidence between the emission line peaks and the infrared 
sources we have identified. Sources of particularly bright 
[S iv] line emission are sources 1 and 4, along with Field 
C. Source 7, on the other hand, is prominent in [Neil] 
emission. The localized nature of the emission allows in- 
vestigations of the physical conditions in the individual 
sources. 

3.3. Electron density 

The ratio of two lines of the same ionization state 
of a single species, emitted from levels with similar 
excitation energies can be used as a tracer of the 
electron density n e (jRubin et al.l I1994I ). We use the 
[SlIl]18.7^m/[SlIl]33.6/xm to investigate the electron 
densities in the NGC 604 region, based on the lores line 
ratios. We use the lores data rather than the hires data 
because the two [S III] lines are measured at similar spa- 
tial resolutions in the lores data, while for the hires data 
they are measured in two different apertures ( [S in] 18.7 



falls on the SH module, while [Sin]33.6 falls in the LH 
module). 

The [S ni]18.7/im/[S ni]33.6/im ratios calculated from 
Table |6] are very uniform in the region, varying from 0.51 
to 1.02, with the highest value near source 7, in the NW 
infrared lobe. These v alues are very close to the low den- 
sity limit discussed in iDudik et all l|2007j ) (their Fig. 9), 
and hence our data provide only approxima te upper lim- 
its for the electron densities. Using the IDudik et al.l 
(2007) diagram corresponding to a gas temperature of 
10 4 K, we derive upper limits for the electron densities 
between log n e = 1.5 cm- 3 and log n e — 2.5 cm- 3 in 
the region. These values ar e in good agreement w i th the 
electron density derived by Maiz-Apclla niz et al.l (2004) 
using the optical line ratio [S n]A6717/A6731. 

3.4. Ionization conditions 

3.4.1. Hardness of the radiation field 

We use the lores [Slv]10.5/xm/[Neli]12.8/im line ratio 
map to visually trace the spatial variations in the hard- 
ness of the radiation field. The use of these two lines has 
an observational motivation. The two lines are within 
the SL wavelength range and hence have the same spatial 
scale and pixel size in the map. This is not the case for 
the [Ne in] /[Neil] ratio, which was a more natural choice. 
The ionization potential required to produce [Ne II] is 
only 21.6 eV, while it takes 34.8 eV to ionize [S in] to 
obtain [Siv]. Hence, the ratio between these two ionic 
species traces the hardness of the radiation field, which 
can be interpreted in terms of stellar ages, with [Neil] 
tracing star formation activity during the last 10 Myr 
and [S iv] tracing massive stars born in the last 4-6 Myr. 
In Fig. [5] we show the corresponding line ratio map ob- 
tained from the spectral cube. We further discuss this 
map in §4.1.31 

3.4.2. Ionization state of the gas 

The mid-IR line ratios do not only depend on 
the hardness of the radiation field. They also de- 
pend on the ionization state of the gas, that can be 
parametrized using the ionization parameter Q (the 
ra tio of the ionizing pho t on de nsity to gas density). 
In iMartmez-Galarza et al.l (|2011l ) (MG11 hereafter) we 
have used the hires li ne ratios as deriv e d fro m the 
measurements done in iLebouteiller et all ((2008), and 
combined them with the radiative transfer models in 
iLevesque et al.l (I2010T) using the interactive ITERA soft- 
ware (|Groves fc Alien! E010), to break the de generacy 
between age and ionization parameter in the partic- 
ular case of the 30 Doradus region. We found in 



TABLE 6 
Lores fine structure lines 3. 



Source 


[Arn]6.9/im 


[Ar ni]8.9/Ltm 


[S iv]10.5/mi 


[Nell]12.8Atm 


[Nem]15.5/im 


[Sm]18.7/mi 


[S lll]33.7Atm 


[Sill]34.8£tm 


Src 1 


0.58(0.19) 


2.05(0.33) 


3.07(0.42) 


4.61(0.64) 


4.83(0.67) 


6.41(0.91) 


6.64(1.13) 


2.54(0.76) 


Src 2 


0.60(0.17) 


0.58(0.16) 


0.84(0.14) 


2.89(0.40) 


2.81(0.40) 


3.53(0.54) 


4.29(0.76) 


1.74(0.53) 


Src 3 


0.59(0.14) 


0.90(0.34) 


1.54(0.28) 


3.18(0.39) 


3.84(0.47) 


5.66(0.83) 


7.90(1.40) 


3.29(0.96) 


Src 4 


0.48(0.16) 


2.05(0.31) 


4.02(0.53) 


3.39(0.52) 


4.08(0.48) 


5.37(0.79) 


6.52(1.13) 


2.57(0.76) 


Src 5 


0.51(0.08) 


0.47(0.09) 


0.32(0.06) 


1.85(0.19) 


1.15(0.13) 


2.21(0.23) 


4.37(0.72) 


1.90(0.47) 


Src 6 


0.15(0.02) 


0.35(0.04) 


0.25(0.03) 


1.28(0.11) 


1.02(0.11) 


1.65(0.01) 


2.34(0.01) 


0.90(0.01) 


Src 7 


0.93(0.21) 


2.59(0.35) 


2.68(0.45) 


6.80(0.86) 


7.02(0.94) 


9.10(1.30) 


8.90(1.11) 


3.15(0.74) 


NGC 604 


0.17(0.04) 


0.29(0.06) 


0.34(0.04) 


1.04(0.11) 


1.04(0.10) 


1.57(0.01) 


1.90(0.01) 


0.77(0.01) 



a Units are 10~ 7 Wm~ 2 sr~ 



TABLE 7 

Fine structure line fluxes as extracted from the hires data. 



Source 


[S iv]10.5/mi 


[Nc Il]12.8^m 


[Nelll]15.5/im 


[SlH]18.7/Ltm 


[S ni]33.7Atm 


[Si n]34.8Atm 




[10~ 20 xW cm" 2 ] 


[10~ 20 xW cm" 2 ] 


[10~ 20 xW cm" 2 ] 


[10~ 20 xW cm" 2 ] 


[10- 20 xW cm" 2 ] 


[10- 20 xW cm" 2 ] 


Field A 


1.63(0.36) 


3.71(0.14) 


4.24(0.08) 


4.54(0.12) 


5.54(0.07) 


1.85(0.09) 


Field B 


2.24(0.22) 


4.24(0.24) 


5.52(0.12) 


5.49(0.12) 


9.03(0.09) 


2.76(0.09) 


Field C 


2.60(0.46) 


5.77(0.24) 


6.63(0.02) 


7.29(0.18) 


10.08(0.09) 


2.71(0.11) 




Fig. 5. — Continuum-subtracted line emission in NGC 604. Color coded are the IRS fluxes in MJy sr 1 . (a) [Siv]10.5/jm emission has 
well defined peaks near sources 1, 4, and field C. (b) [Nell] peaks near source 7. The white contours are IRAC 8 /an emission, and the 
RMS noise values are 0.45 MJy sr -1 and 4.5 MJy sr -1 , respectively. 



that study that the [Ne iii]15.5/um/[Neri]12.8/zm and the 
[Srv]10.5/xm/[Siii]18.7/zm ratios in 30 Doradus are com- 
patible with starburst ages < 3.0 Myr and ionization pa- 
rameters 8.0 cm s _1 < log Q < 8.6 cm s . 

We perform a similar analysis for NGC 604 using the 
lores lines, which are extracted from the same spatial 
area regardless of the IRS module in which they fall. 
In Fig. we plot line r atios computed from Table H] 
as compared to a set of iLevesaue et ail (|2010l ) models 
of an instantaneous starburst with sub-solar metallicity 
(Z = QAZq) and low electron density (n e = 10 cm -3 ), 
in accordance with the low density regime inferred from 
the line ratios. The overall line ratios are indicative of 
ages between 4 and 4.5 Myr for the individual sources. 
Although both the observational uncertainties in the line 



ratios and the spatial confusion between the unresolved 
sources in the IRS map affect the measured line ratios, 
Fig. indicates a similar age for all the infrared bright 
sources. However, sources 1 and 4 show larger ionization 
parameters, sources 2, 3, and 7 show moderated ioniza- 
tion parameters, and sources 5 and 6 have the lowest 
ionization parameters. While this is not necessarily in- 
dicative of a stronger radiation field near sources 1 and 
4, it points to a large number of ionizing photons near 
these sources. 

Using the hires line fluxes we also measure the strength 
of the radiation fi eld in fields A, B, and C with the 
parametrization of Bci rao et al.l (|2006f) . which uses the 
[Ne in] /[Nell] to estimate the strength as the product of 
the field intensity and the field hardness: 
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Fig. 6.— The [Slv]10.5/mi/[Nell]12.8/im ratio in NGC 604 as 
derived from the spectral map. Color coded is the value of the ratio. 
The RMS noise equals 0.12. For reference, the IRAC 8 fim contours 
are superimposed as the black solid line. The measured ratios are 
used to estimate the hardness of the radiation field towards specific 
locations (see text). 




1.5 -1.0 -0.5 

log„( SIV 10.522um/ SHI 18.682um) 

Fig. 7. — Compari s( 
ILevesque etaH j20ToT) models. The [Nem]15.5/mi/[Nen]12.8/im 
ratio is plotted versus the [S lv]10.5/im/[S lll]18.7/im ratio. Data 
points from Table[6]are overplotted to the grid of models, that cover 
a range of ages from 4 Myr (blue) to 4.5 Myr (cyan) and a range of 
ionization parameters from 8 X 10 7 cms -1 (red) to 4 X 10 8 cms -1 
(yellow) . The line ratio for the integrated map is indicated by "N" , 
to the upper right of the symbol. Based on the propagation of the 
flux uncertainties listed in Table [6] typical errors in the shown line 
ratios are 0.1 dex. 



[NcII]12.8/jm 



F 



F 



[NcIII]15.6/imJ 
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We get a radiation field strength of (9.1 ± 0.5) x 
10- 20 W cm- 2 , (12.7±0.9) x 10- 20 W cm~ 2 and (14.25± 
0.5) x 10~ 20 W cm -2 for fields A, B, and C, respectively. 

3.5. [Sill] emission 

Most of the lines listed in Table [6] are from regions with 
ionized hydrogen gas. The exception is [Si II], which orig- 
inates in a variety of environments including H II regions , 
but also X-ray dominated regions ()Malonev et al.lll996l) . 
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Fig. 8. — The [Sill]/[Nell] ratio plotted against the ratio of 
PAHi7 Mm strength to the strength of all the other PAH features, 
for our sources. Error bars in both quantities are shown. 



TABLE 8 

H2 ROTATIONAL LINES 



Source 


H 2 S(2)12.27Atm 


H 2 S(l)17.03|tm 




[10" 20 xW cm~ 2 ] 


[10- 2O xW cm" 2 ] 


Field A 


0.19(0.15) 


0.20(0.01) 


Field B 


0.13(0.15) 


0.13(0.02) 


Field C 


0.16(0.04) 


0.14(0.01) 



high density PDRs (Kaufman ct al. 2006), and regions of 
shocked gas, where heavy elements are returned to the 
gas phase. It is in general hard to pin down the physical 
mechanism for the emission of [Si 11] . We investigate this 
in NGC 604 by looking at a possible correlation between 
the ratio [Si 11] /[Neil] and the ratio of PAH emission at 
17 /xm to all the other PAH features together. We plot 
the relation between these two ratios in Fig. [8] In £ 14.2.71 
we discuss this finding in the context of several possible 
scenarios for the emission of [Si 11] . 

3.6. H2 emission 

Additional diagnostics on the physical conditions in 
NGC 604 come from the molecular hydrogen lines aris- 
ing from pure rotational transitions at 12.27 /zm (0-0 
S(2)) and 17.03 (0-0 S(l)), which we detect (although 
marginally in the case of the 12.27 /im line) in the hires 
modules of the IRS in fields A, B, and C. Using the same 
procedure as for the nebular lines, we fit Gaussian profiles 
to the H 2 rotational lines and estimate their strengths. 
In Table [5] we list the measured line strengths. 

H2 temperatures can be calculated from the ratio of 
the two detected line strengths listed in Table [8j us- 
ing a general meth od (see, for example the Appendix of 
iBrandl et aIll2009D . This method holds under the follow- 
ing assumptions: (i) the gas is in local thermodynamic 
equilibrium (LTE); (ii) the hydrogen rotational lines are 
optically thin; (in) the critical densities for the lines are 
< 10 3 cm -3 ; and (iv) the two states of the H2 molecule, 
ortho-H2 and para-H2, exist in a ratio of 3 to 1. From 
the Boltzmann statistics that describe the distribution 
of energy states for the H2 molecule and the relation be- 
tween the line strength and the H2 column density, it can 
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be shown that the excitation temperature is given by: 

_ E 2 -E 1 

where El and E2 are the energies of the involved lev- 
els, k is the Boltzmann constant, C is a constant depen- 
dent on the ratio of Einstein coefficients and statistical 
weights, Ai and A2 are the rest-frame wavelengths of the 
observed lines, and F\ and F% their respective integrated 
line fluxes. Fields A and B have uncertainties close to 
100% in the H 2 S(2) flux (Table ©, and hence we are 
unable to say anything conclusive about the gas tem- 
perature. Field C is comparatively better, and for this 
source we obtain T ex between 200 K and 800 K. 

3.7. SED modeling 

In MG11 we have presented a Bayesian fitting tool to 
study the infrared SEDs of star-forming systems. The 
tool is based on a grid of the Dopita & Groves mod- 
els (D&G models hereafter) , which are thoroughly d e- 
scribed in a series of papers (iDopita et alJl2005l [20061 13: 
iGroves et al. 2008). Here we briefly describe the tool and 
then apply it to the observed infrared SED of NGC 604 
in M4.21 including both the IRS spectrum and the PACS 
photometry, to derive statistically meaningful values for 
the following physical parameters of the region: cluster 
age and stellar mass, compactness, fraction of total mass 
contained in embedded objects and fraction of the lumi- 
nosity arising in photon-dominated regions. 

3.7.1. Physics 

The D&G models combine stellar synthesis, radiative 
transfer, dust physics, and self-consistent dynamical evo- 
lution of individual Hn regions to simulate the UV to 
sub-millimeter SEDs of star-forming regions, including 
thermal emission from dust and PAHs, and nebular line 
emission. Using all the available spectral information, for 
a given metallicity (Z) and ambient interstellar pressure 
(Po/k), the fitting tool computes probability distribu- 
tion functions (PDFs) for key model parameters such as 
stellar cluster age (t c \), cluster mass (M c \), compactness 
(C), fraction of the total luminosity arising in the PDRs 
(/pdr), and the mass contribution from a component of 
young embedded objects that we model as Ultra Com- 
pact Hn Regions (M em b). 

We have defined these parameters in MG11. Of rele- 
vance here is the d efinition of C, originally described in 
IGroves et al.l (|2008h . This parameter controls the form of 
the FIR continuum and the wavelength of the FIR peak. 
The shape of the SED at wavelengths longer than about 
15 /jm is a function of the distribution of dust temper- 
atures throughout a given region. Such distribution is 
a function of the total heating flux incident on the dust 
grains. For a spherical expanding nebula, the heating 
flux is calculated as L^/R^ ill at each time, where L* is 
the luminosity of the cluster and i?nn * s ^ ne ra( hus of the 
Hn region. Hence, models that preserve the run of dust 
temperature with time should also preserve the quantity 
< L*(t) > I < -Rffli(i) 2 >• The compactness parameter 
is proportional to this ratio, and hence sets the evolu- 
tion dust temperature with time. Intuitively, for a given 
cluster luminosity, denser H 11 regions will have smaller 
radii, and hence hotter dust. 



3.7.2. Priors 

Bayesian inference allows to include any previous evi- 
dence on the values of the model parameters in the cal- 
culation of the posterior PDF. Here we use flat, bounded 
priors for all the model parameters, with boundaries set 
by observational and theoretical studies of star-forming 
regions, as discussed in MG11. We assume solar metal- 
licity and adopt a thermal pressure of the surrounding 
interstellar medium (ISM) of Po/k = 10 5 K cm~ 3 . 

The assumption on ISM pressure is supported by mea- 
surements of far-IR line ratios fro m the ISO satellite o n 
a sample of star- forming galaxies (jMalhotra et alJfeOOlh . 
The metallicity of NGC 604 has be en measured to be 
about half-solar (jMagrini et alJl2~007j ). However, our ex- 
periments show that only the solar metallicity models in 
our grid are able to reproduce the observed ratio of PAH 
strength to continuum emission at 100 in NGC 604. 
The discrepancy is most likely due to the specific PAH 
template used in the D&G models, which is an empiri- 
cal template based on observations of starburst galaxies. 
The choice of metallicity does not only affect the PAH 
emission strength. At far-IR wavelengths, the change in 
dust column and mechanical luminosity of the starburst 
with metallicity produces a slight shift and a broadening 
of the far-IR bump. While this implies an additional de- 
generacy with the compactness parameter, the additional 
errors introduced in the determination of C are smaller 
than our chosen step size for the compactness grid, which 
is 0.5. 

3.7.3. x 2 Weighting 

In MG11 we have included a discussion of the observa- 
tional uncertainties involved in measuring the IRS fluxes, 
which include absolute and relative flux calibrations and 
systematic errors due to specific observing conditions. 
Based on that discussion, for fitting purposes we have 
adopted a uniform uncertainty of 10% across the IRS 
wavelengths, which also implies a uniform weighting for 
all data points in the fitting. 

The x 2 minimization procedure described in MG11 has 
been slightly modified here to ensure that the x 2 mini- 
mization is not dominated by the IRS range, which has 
many more resolution elements, as compared to only two 
data points at 100 /im and 160 /im from PACS. Each bin 
contributes to the % 2 with a weight that is proportional 
to the bin size in the logarithmic wavelength space. The 
bin size is set by the resolution of the models in the IRS 
range and by the wavelength separation between data 
points for the PACS data. This results in a weighting 
function that increases uniformly in the IRS range so 
that the bins at the short wavelength end (around 5 fim) 
have about half the weight of those at 35 /im and about 
0.25 times the weight of the PACS bins. In section S R2I 
we apply this tool to the observed SED of NGC 604. 

3.7.4. Color Correction of PACS photometry 

The D&G models compute the monochromatic flux 
densities for each wavelength bin. However, the flux 
densities measured by the PACS filters at 100 ^m and 
160 fim are not monochromatic. They are the integrated 
flux densities over certain wavelength ranges, as modu- 
lated by the filter response function. Before we compare 
the model fluxes to the observed PACS photometry, it is 
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necessary to evaluate the errors introduced by this dif- 
ference. To do so, we adopt the method described in 
Ida Cunha et al.l (|2008|) to predict the flux density of a 
source when observed using a given filter. According to 
their approach, the flux density for any filter is: 



F 



\2 



J dA F x XR X 
J dX C A A R a 



where 



An = 



/ dA A R A 
J dA R A 



(3) 



(4) 



is the effective wavelength of the filter response R a , c is 
the speed of light, and C\ is a calibration spectrum that 
depends on the photometric system used for the calibra- 
tion. In the case of the PACS photometer, the calibration 
spectrum is constant (CaA = constant), which simplifies 
the calculations. Using this method, and the filter re- 
sponse curves available at the Herschel Science Center, 
we have computed effective wavelengths for the PACS 
filters at i* 00 = 102.6 fim and L* 60 = 167.2 fim. Fur- 
thermore, using Eq. [3j we have estimated by how much 
the PACS fluxes, as predicted by our models, should be 
corrected once the color correction is applied. For our 
best fit model, the differences between the monochro- 
matic model fluxes and those that would be measured 
by PACS are < 1% for the PACS green filter (100 H 
and - 9% for the PACS red filter (160 /im). These un- 
certainties are within the 10% observational errors. 

4. DISCUSSION 

4.1. Notable Sources 

The most striking morphological characteristic of the 
sources labelled 1-7 in Fig.[3]is that most of them (sources 
1-5) are well defined, individual infrared-bright knots. 
Sources 1, 2, 7, and 4 are, in that order, the strongest 
sites of PAH emission, and also the sites with warmer 
dust continuum, as shown in Tables [3] and [4] This implies 
warmer dust temperatures near sources with stronger 
8 /im emission, where most of the MYSO candidates 
have been identified by iFarina et al.l (|2012D using NIR 
photometry. 

4.1.1. Source 2 

This source, located in the SE lobe, is bright in all 
IRAC bands but has no optical, Ha or FUV counter- 
part. In addition, it has one of the strongest silicate 
absorptions at 10 /nn as shown in Fig. 02 and it is only 
about 5" fr om the peak of o n e of t he CO clouds first 
reported in IWilson fc Scovillel (jl992ft . Using HST data, 
iMai'z-Apellaniz et al.l (|2004H derived an extinction map 
for the region. They found a strong peak of extinction 
at the location of source 2, which is consistent with its 
relatively strong silicate absorption. Compared with all 
the other sources, nebular emission towards this source 
is weak (Table Its high optical extinction, bright 
PAH emission, and spatial coincidence with a reservoir 
of molecular gas makes of source 2 a very good candidate 
for a site of embedded star formation. 

4.1.2. Source 5 



Source 5 shows the coldest dust temperature among 
our sources, as traced by the Fi 5lirn /F 30lim ratios 
listed in Tabled Its [Siv]10.4/mi/[Siii]18.7/im and 
[Nein]15.5/im/[Neii]12.8/im line ratios derived from Ta- 
ble|6]are both the lowest among our sources and the PAH 
features listed in Table [4] are weaker towards this source, 
compared with all the other well defined IR knots. The 
combination of cold dust, low ionization state, soft radia- 
tion field (as traced by the line ratios) , and weak emission 
from PDRs are indicative of a more evolved stage for this 
particular source. 

Source 5 is also brighter than all the other sources 
at 3.6 /im. This band traces the 3.3 /im PAH feature, 
and photospheric emission from stellar populations older 
than 10 Myr. Since the continuum PAH emission is weak 
towards this source, the most plausible explanation for 
the excess at 3.6 /im is again the presence of stars older 
than 10 Myr. In fact, this source coincide s with the lo- 
cation of one of the WR stars identified bv lHunter et al.l 
(1996), and is located near the boundary between the ac- 
tive star-forming western and quiescent and older eastern 
hemis phere of NGC 604, as described bv lTullmann et all 
(l2008h . 

4.1.3. Source 7, Source 4, md Source 1 

Source 7, located in the NW infrared lobe, shows the 
strongest nebular line emission among our sources (Table 
[5]). This source is very close to a bright ridge of photo- 
spheric optical emission from a nearby group of young 
stars (see Fig. [2]). The higher electron density near the 
NW lobe is consistent with the strong nebular lines mea- 
sured near source 7 and with the enhanced Ha emission 
observe d in this same area of NGC 604 (see, for example, 
Fig. 2 in lRelano fc Kennicuttll2009| ) . and implies a higher 
ionization state. This source is most likely the location 
of the youngest main sequence stars of the cluster. This 
is supported by our measurements of the radiation field 
strength in the neighboring field C, where we have mea- 
sured a stronger radiation field using Eq. [Q relative to 
the other fields. In fact, UV spectroscopy using HST's 
Space Telescope Imaging Spectrograph has revealed that 
field C coincides with the location of a young, luminous 
07 star, most likely one of the most massive stars in 
NGC 604 (|Bruhweiler et al.ll2003] ). 

The [Siv]/[Neli] ratio map of Fig. [6] peaks within the 
projected area of source 4, which has a relatively harder 
radiation field, where we have also measured a warm 
dust component from the spectral continuum slope (Ta- 
ble [3]). Another locat ion with relatively hard ra diation 
field is source 1. The iRelano fc Kennicuttl (|2009ft study 
of the region reveals enhanced Ha emission in sources 
1, 4, and 7. However, the eastern part of NGC 604, 
where sources 1 and 4 are located, has a higher op- 
tical extinction traced by the Balmer optical thickness 
(■?Bai ~ 1-2), as compared w ith the surrounding averag e 
extinction (r Ba i = 0.2-0.3) (iMaiz-Apellaniz et al.ll2004|i . 
The [Siv]10.5/xm falls on the broad silicate absorption 
feature near 10 /im, and this leads to an underestimation 
of the [S iv] / [Ne n] ratio in regions of high extinction. Al- 
though this does not affect our general picture about the 
radiation hardness. Larger [Siv]/[Neii] ratios and hence 
younger ages might be expected in sources 1 and 4 from 
extinction corrected [Siv]/[Neii] ratios. 

Our results are consistent with the optical spectroscopy 
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Fig. 9.— Best fit model to the integrated SED of NGC 604. The 
IRS and PACS data are shown as red diamonds, and the IRAC 
photometry as blue diamonds. The solid line is the best fit model. 
Also shown are the contributions from pure H II region emission 
(dashed blue), PDR (solid blue), embedded population (dotted 
blue), and photospheric emission from stars older than 10 Myr 
(dashed red). Residuals are shown in the bottom panel. 

of the region presented in iMai'z-Apellaniz et al.l (|2004f ) 
in terms of the high density and excitation derived near 
sources 1, 4, and 7, as well as field C, and point to the 
fact that these are the regions where the youngest, most 
massive stars are located within NGC 604. Finally, the 
[Slv]10.5/mi/[Neii]12.8/xm ratio is sensitive not only to 
the cluster age, but also to the gas density. However, 
at the relatively low electron densities measured in the 
region (logn e ~ 1.5 — 2.5 cm -3 ), we do not expect the 
measured line ratios to be significantly affected. 

4.2. The evolutionary status of NGC 604 

Fig. [9] shows the best fit to the observed integrated 
SED of NGC 604 (IRS + PACS) using the Bayesian 
tool described in £13.71 Also, in Fig. [10] we show the 
resulting PDFs, covering the totality of our parameter 
space. In Table [9] we list the associated best fit param- 
eters and uncertainties, with the values for t c \ a n d M c \ 
determined independently bv lEldridge fc Relanol (|2011f ) 
using several methods, including optical SED fitting. We 
also include / em b, which is the ratio of mass contained 
in young embedded objects (M em b) to stellar mass in 
the cluster (M c \). This Bayesian fit, as well as the other 
spectral features discussed above, provide a wealth of in- 
formation about the physics and the evolutionary stage 
of NGC 604. 

4.2.1. Dust temperature and compactness 

The thermal continuum, PAH emission, and line emis- 
sion are well reproduced by the best fit. The com- 
bined IRS and PACS data are consistent with a far-IR 
SED peaking at approximately 70 /im, which indicates 
an effective dust temperature of ~ 40 K. Using multi- 
wavelength photom etry of a sample of H n regions in the 
Magellanic Clouds, lLawton et al.l (|2010f) show that the 
infrared SEDs of most of these regions peak at 70 /im. 
Our result indicates that this average dust temperature 
is not unique of "local" Hn regions, but also applies to 
at least one relatively more distant region. 
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a Parameters listed are cluster age (t c i), compactness (logC), 
fraction of total luminosity arising in PDR regions (/pdr)> 
stellar mass of the cluster (M c \), mass contained in embedded 
objects (M om b), and ratio of embedded to stellar mass (/ em b)- 

We have stated that the compactness parameter con- 
trols the incident heating flux L^/Rhii as a function of 
age of the Hn region. From our given set of best fit 
parameters it is then p ossible to derive an incident heat- 
ing flux (see Fig. 3 in iGroves et al.l I2008T ). We obtain 
an incident flux of log(L*/i?Hii) = 0.65 erg s _1 cm~ 2 for 
our best fit values of compactness and age. The bolo- 
metric luminosity of the cluster can be obtained from 
integration of the best fit SED, and gives approximately 
6.0 x 10 41 erg s _1 , which allows us to solve for the effec- 
tive radius of the Hn region. We obtain i?nn ~ 120 pc. 
This value corresponds to the scale of the larger filaments 
observed in NGC 604 (Fig. [2]). This agreement between 
a model-derived quantity and a measurable observable 
demonstrates the capabilities of our Bayesian approach 
to obtain valuable physical information when applied to 
unresolved star forming regions. 

4.2.2. Stellar mass 

The total stellar mass of l-6t}'g x 1O 5 M that we esti- 
mate with the SED fitting is in agree ment, within the un- 
certai nties, with the mass derived by Eldridge fc Relanol 
(2011). Besides, within the uncertainties, this mass is 
similar to the mass derived for 30 Doradus in MG11, 
which equals 0.63±£;?J x 10 5 M Q . 

In general, previous studies have found that the total 
stellar mass in 30 Doradus is larger than in NGC 604. We 
propose two reasons to explain the fact that the value of 
our derived mass for NGC 604 is larger than the value 
of the mass derived for 30 Doradus. First, the physical 
projected area covered in our NGC 604 map is relatively 
larger than the projected area of the 30 Doradus map, 
and hence there is some missing emission from the pe- 
riphery in the latter case t hat do es not enter the SED. 
Second, we observe in Fig. 10(c) a degeneracy between 
the cluster mass and the fraction of the total luminosity 
from PDRs that covers a range of masses of ~ 0.5 dcx. A 
similar degeneracy was found for 30 Doradus in MG11. 
We have argued that the reason for this degeneracy is the 
fact that both the PDR and the pure H II region spectra 
contribute to the dust thermal continuum. An increase in 
PAH from emission is accompanied by an increase in the 
overall infrared luminosity. This also renders the addi- 
tional far-infrared data insufficient to solve this degener- 
acy. Optical wavelengths data, where the emission from 
PDR is significantly different from that of an H II region, 
are needed to settle this issue. Nevertheless, although our 
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result is subject to this M c \ — /pdr degeneracy, a compar- 
ison between our derived masse s and the masses de rived 
for 30 Doradus and NGC 604 in lHunter etUI (|1995l ) and 
lEldridge fc Relanol (|2011h respectively, implies that our 
Bayesian fitting tool can find solutions that are in rea- 
sonable agreement with published literature values. 

A brief discussion o n the mass of the individual in- 
frared sources in Fig. 2(b) is relevant here. Although 



Fig. 10. — Two-dimensional PDFs for selected pairs of parameters covering the totality of the parameter space. The pairs shown are: (a) 
logC - t c i, (b) logM cm b - logM c i, and (c) /pdr, - logM c i. The color code indicates the normalized probability. The cross symbols mark 
the best-fit values while the color contours indicate the 1-cr (blue) and 90% (green) confidence levels. 

at longer wavelengths, diffuse dust heated by older stellar 
populations starts to dominate the SED. In the following 
we discuss how our best fit parameters for NGC 604 help 
us constrain the age of its dominating stellar population. 

The fraction of total mass contained in embedded ob- 
jects that we measure here for NGC 604 (/ cm b ~ 0.08) 
is smaller than that measured for 30 Doradus in MG11 
(/cmb ~ 0.53). If we interpret this result in terms of re- 
cent star formation taking place in these two regions, this 
suggests that 30 Doradus has a larger relative amount of 
mass in embedded objects, but this conclusion should be 
taken with caution, since other dust heating mechanisms 
can play a role. Additionally, the age of ~ 4 ± 1 Myr de- 
rived with our fitting method is consistent with the study 
of the line ratios that w e presented in 33.41 and a grees 
well with the studies by lEldridge fc Relanol ()2011l ) and 
iHunter et all (|l996l) . In MG11 using the same method, 
for 30 Doradus we have obtained an age of 3 ± 1.5 Myr. 
Our results suggest that the dominating stellar popula- 
tion in each of these two giant HII regions (i.e., the one 
that provides most of the ionizing photons) are consis- 
tent with slightly different model ages when variations in 
all model parameters are considered 



the SED models are not intended to reproduce ensem- 
bles of embedded objects but entire Hll regions, we fit 
their spectra using our tool to estimate their individual 
masses. In Fig. [TT] we show the fit to the IRS spectra 
of two representative sources that significantly differ in 
their IRAC colors: sources 1 and 5. Note that the masses 
of these sources are of the order of a few times 10 3 Mq. 
This is the case for all the other sources, which are sim- 
ilarly luminous. Although we do not attempt a detailed 
interpretation of the fitted parameters, because we are 
aware of the limits of our models, it is qualitatively in- 
teresting that we derive a higher PDR contribution and 
a higher fraction of embedded mass for source 1 than for 
source 5, which is consistent with the latter being in a 
region dominated by stars older than 10 Myr (see §4.1|) . 

4.2.3. Age of the region 

Both the instantaneous burst approach that leads to a 
single, coeval stellar population and the continuum star 
formation rate approach that leads to a uniform distri- 
bution of stellar ages, are opposite idealizations of the 
star formation history in a particular starburst region. 
In fact, rather than showing one of these limiting situa- 
tions, it is often found that several star formation events 
in a particular region lead to different stellar popula- 
tions of d ifferent ages. NGC 604 is not an exception to 
this rule. lEldridge fe Rclanc] (|2011h argue that the pop- 
ulations of Wolf-Rayet (WR) stars and red supergiants 
(RSG) in NGC 604 belong to different formation events, 
with estimated ages of 3.2 ±1.0 Myr and 12.4 ±2.1 Myr, 
respectively. A similar situation is observed in other gi- 
ant Hn regions such as 30 Doradus. 

Henceforth, when we refer here to the age of the region 
what we mean is the age of the dominating stellar popu- 
lation, which is usually composed of the youngest, most 
massive main sequence stars in these systems. As we 
have already discussed, most of the infrared continuum 
and line emission up to wavelengths of 100 /im, arises 
from the interaction of the strong radiation field from 
young (< 10 Myr) stars and the surrounding ISM. Only 



4.2.4. PDR content 

Although the fraction of total luminosity arising in 
PDRs is dege nerate with the total cluster mass, as shown 
in Fig. 10(c)| our result from SED fitting suggests that at 
least half of the energy output from NGC 604 is gener- 
ated in PDRs, the other half arising in ionized regions 
closer to the cluster stars. The presence of a lumi- 
nous PDR in NGC 604 is consistent with the findings of 
iHeiner et al.1 ((2009) , who find considerably high amounts 
of atomic hydrogen (Njjj ~ 2.0 x 10 21 cm -2 ), and asso- 
ciate them with the photodissociation of H2 in a dense 
(n ~ 500 cm -3 ) PDR. 

4.2.5. Molecular hydrogen in NGC 604 

The H2 excitation temperature near field C (T cx ~ 
200 K — 800 K) , calculated using Eq. [5] is indicative of 
the average temperature of the warm gas in NGC 604. 
Additionally, for a given value of the molecule angular 
momentum, the total number of hydrogen molecules is 
proportional to the strength of the line. 

Our results on the H2 S(l) line, which has smaller 
uncertainties (Table [8]), indicate that warm molecular 
hydrogen is present in fields A, B, and C and slightly 
more abundant towards field A, which has the strongest 
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Fig. 11. — SED fits to the spectra of source 1 (a) and 5 (b). Color 
consistent with masses of a few times 10 3 Mq. 

H2 lines. This res ult is consistent with the findings of 
iTosaki et al.l ([2007), who report the detection of an arc- 
like distribution of warm and dense molecular gas that 
extends across coincident our fields A, B and C, and as- 
sociate it with gas compression by the expanding H II 
region and sequential star formation in NGC 604. The 
center of Field A is, in projection, about 5" (or 20 pc) 
away from the peak of CO emission first reported in 
IWilson fc Scovillel (11992D and c o nfirme d to be a dense 
molecular cloud in iMiura et al.l (j2010t ). although given 
the extension of this cloud (about 30 pc), we can consider 
it as coincident with Field A. We cannot confirm a phys- 
ical association of the molecular cloud with our infrared 
sources, since we only see a projected spatial association. 
However, based on our results and the mentioned previ- 
ous work, it becomes evident that both dense, molecular 
gas as well as warm molecular hydrogen is detected in 
the line of sight towards Field A. 

4.2.6. Selective destruction of small PAHs 

The 17 fim PAH feature is generally associated with 
out-of-band bending mode s of large PAH molecules, con - 
taining w 2000 C-atoms (jVan Kerckhoven et all 12000( 1 . 
Table [5] shows that the relative strength of this feature 
exhibits an increasing trend with the X-ray flux. If we 
adopt the PAH size argument, this implies that it is the 
ratio of large to small PAH molecules that scales with the 
X-ray field. In fact, Fig. shows that source 3, which 
shows the largest enhancement of the 17 /xm PAH feature 
among our sources, coincides in projection with a peak 
of soft X-ray emission. Soft X-ray emission from mas- 
sive stars can be associated with shocked stellar winds 
or magnetically co nfined gas near the wind base, near 
the stellar coronae (jCassinelli fc Swan kl ll~983l l. 

The disso ciation of PAH molecu les by X-rays has been 
discussed in lMicelotta et al.l (|2010H . They argue that not 
all X-ray photon absorptions by PAH molecules lead to 
photodissociation and estimate that after a second elec- 
tron has been ejected by the PAH molecule via the Auger 
effect, the molecule is left with an internal energy of 14- 
35 eV, enough to dissociate small (50 C-atoms) PAHs, 
but possibly insufficient to dissociate larger molecules. 
They leave the question of large PAH survival open. Our 
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Fig. 12.— Soft X-ray emission in NGC 604. The black and white 
image is the IRAC 8 fira emission. The red contours correspond to 
the ACIS image. Source 3 coincides with a local peak of soft X-ray 
emission with a flux of 3.70 X 10 — 9 erg cm -2 s — 1 (see Tablc[5]|. 



observations suggest that, at least in the particular case 
of source 3, X-ray emission from massive stars may be 
responsible for the dissociation of small PAH molecules, 
creating an enhancement of the 17 /jm feature. Similar 
relative enhancements of the 17 /im feature have been ob- 
served in regions of har d radiation fields in galaxies, such 
as th e vicinity of AGNs ((Smith et al.ll2007HO'Dowd et all 
2009). It is worth mentioni ng here that a MY SO candi- 
date has been identified by (iFarina et alll2012ri le ss than 
1" away from source 3. Also. iBarba et al . (2009) report 
the presence of at least 2 red supergiants (RSGs) coinci- 
dent with the location of this source. If other RSGs have 
already exploded as supernovae near this location, this 
could provide part of the X-ray enhancement observed. 

4.2.7. The origin of the [Sill] emission 

Fig. [8] is suggestive of a correlation between the en- 
hancement of [Sin] and the enhancement of 17 /im PAH 
emission at scales of tens of parsecs. Despite the rela- 
tively large uncertainties in both the PAH ratio and the 
[Si n]/[Ne n] ratio, the figure indicates clearly that the 
enhancement of the PAH 17 /im feature scales with the 
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relative strength of the [Si n] emission, with the exception 
of source 5, which lies outside this correlation. This re- 
sult suggests a common underlying physical mechanism 
for the enhancement of both features. We have shown in 
Table [5] that the regions of strongest 17 /im PAH emission 
also correspond to regions of enhanced X-ray emission. 
Based on the triple correlation between PAH^m emis- 
sion, [Sin] emission, and soft X-ray flux, we argue that 
absorption of ambient X-rays may be linked to both the 
destruction of small PAH molecules and the excitation of 
[Si n] , and produces the correlation observed in Fig. HI 

As we have pointed out, source 5 is an outlier of this 
correlation: it shows a relatively large [Sin] /[Neil] ra- 
tio, but very moderate PAHi7 Mm /EPAff ratio. It also 
lies in a region of low X-ray flux, near the edge between 
the western, active star-forming hemisphere of NGC 604 
and the quiescent eastern region where most o f the Wolf- 
Rayet stars are located (|Drissen et al.l 12008). In fact 
it coincides with th e location of one of these WR stars 
([Hunter et all 11996( 1 . In this part of NGC 604, the rel- 
atively high X-ray luminosity is most likely powered by 
the shocked gas resulting from older dynamical processes. 
Based on these facts, we speculate that the [Siii]/[Neii] 
ratio for this source is enhanced due to the presence of 
shocked gas rather than X-ray induced. Although no 
supernova remnants (SNRs) have been identified in the 
eastern hemisphere, this is not surprising, since such ob- 
j ects expanding into a l ow-density gas are hard to detect 
dChu fc Mac LowlllMih . 

4.3. Ongoing and Triggered Star Formation in 
NGC 604 

Several authors have gathered evidence that sug- 
gest the presence of ongoing star formation in 
NGC 604 (iMai'z-Apellaniz et alj|2005 iTosaki et al.| [2007l: 
iRelano fc Kennicuttj|2009HFarina et al.ll2012H . Neverthe- 
less, none of these authors have provided spectroscopic 
confirmation for such very recent star formation activ- 
ity in the region. Our study of the sub-condensations 
whose spectra are shown in Fig. [3] allows us to quantify 
the amount of embedded star formation happening in 
NGC 604. 

In 34.11 we have discussed the IR-bright sub- 
condensations that we detect in the IRAC maps. At 
the distance of NGC 604 it is unlikely that they are in- 
dividual massive young stellar objects. Their SED lumi- 
nosities are consistent with stellar masses of a few times 
10 3 Mq (Fig. [II]). The exact values depend on the indi- 
vidual SED fits and the uncert ainty in the determination 
of the stellar mass (Fig. 10(b) ). This implies that all to- 
gether they account for between 5% and 20% percent 
of the total stellar mass in NGC 604. We have shown in 
Fig.[7]that the line ratios measured towards these clumps, 
are consistent with ages of between 4 Myr and 4.5 Myr. 
Although the timescales for the dissipation of dusty en- 
velopes around star-forming clusters can be as long as 
10 Myr in regions of high density, the low electron den- 
sity that we have inferred in §3.31 implies that NGC 604 
is a rather diffuse region. Consequently, it is unlikely 
that our infrared clumps are evolved clusters. The pres- 
ence of these young, massive, IR-bright clusters, some of 
which have deep silicate absorption features, is a clear 
indication of significant massive star formation currently 
taking place in NGC 604. 



Additional evidence that supports this scenario comes 
from our SED fitting analysis. Fig. [9] shows that a warm 
component of dust is necessary to fit the spectrum of 
NGC 604 between 15 //m and 30 /im. This component 
(dotted blue line in Fig. [9]) arises from MYSO (that in 
the D&G frame we model as Ultra Compact Hn Re- 
gions), and corresponds to dust temperatures of about 
300 K. Although the required fraction of mass in em- 
bedded objects is considerably smaller than that derived 
for 30 Doradus in MG11, it suggests the presence of a 
considerable amount of embedded star formation in the 
region. The total embedded mass derived is in fact com- 
parable to the sum up of the individual clump masses, 
and equals approximately 8% of the total cluster mass, 
as can be calculated from Table [9j 

In an evolutionary context, sequential star formation 
in NGC 604 is a plausible scenario to explain our obser- 
vations. This idea has been explored by other authors 
before. Based on sub-millimeter observations o f the CO 
(J = 3-2)/CO (J = 1-0) ratio. ITosaki etHl ([2001 re- 
port the existence of a dense ridge of molecular gas that 
surrounds the main cluster in NGC 604 and extends in 
the SE-NW direction, closely following the location of 
our bright IR lobes. To explain their results, they adopt 
an scenario in which the compression of molecular gas 
by the mechanical input from the main cluster (the first 
generation of stars) has triggered a second generation 
of stars near the NW infrared lobe. The strong radia- 
tion field that we observe close to field C and source 7, 
where main sequence stars are clearly observed, supports 
the existence of this second generation of highly ionizing 
stars. Furthermore, the results discussed in this sub- 
section indicate that massive star formation is currently 
taking place within the SE lobe of NGC 604, even further 
away from the main cluster, where we have identified the 
IR-bright sub-condensations of Fig. [2j some of which are 
heavily enshrouded by dust (e.g., source 2). 

5. CONCLUSIONS 

We have investigated the physical conditions and quan- 
tified the amount of ongoing massive star formation in 
the star forming region NGC 604. We used a combina- 
tion of observational and modeling tools, including in- 
frared spectrophotometry and Bayesian SED fitting of 
Spitzer and Herschel data. Here are our main findings: 

1. We have identified several individual bright in- 
frared sources along the luminous PDRs that sur- 
round the ionized gas in NGC 604. These sources 
are about 15 pc in diameter and have luminosity 
weighted masses between 10 3 Mq and 10 4 M . 

2. The deep 10 /im silicate absorption feature, mid- 
IR continuum slope, and atomic line ratios towards 
some of these sources indicate that they are young 
embedded systems, and most likely the sites of on- 
going massive star formation in NGC 604. Some 
of them (i.e. source 2) are also associated with gas 
reservoirs as traced by CO maps. This is in agree- 
ment with previous studies of the region. 

3. A scenario of ongoing massive star formation in 
NGC 604 is supported by Bayesian fitting of the in- 
tegrated spectrum (lines+continuum) of the region 
constructed from Spitzer-IRS and Herschel-PACS 



16 



observations. Our results indicate that embedded 
star formation can account for up to 8% of the total 
stellar mass in NGC 604. 



4. The spectral fitting also implies an age of 4.0 ± 
1.0 Myr for the region and a total stellar mass of 
~ l-6_ilo x 10 5 Mq. These results are in agreement 
with independent measurements of these quanti- 
ties using optical broad band photometry. Addi- 
tionally, our best fit model implies an average dust 
temperature of ~ 40 K. 

5. We measure a stronger than average radiation field 
near our sources 7 and field C. This result is con- 
sistent wit h the sequent i al sta r formation scenario 
adopted in lTosaki et al.l (|2007l ). according to which 
a second-generation of star formation in the loca- 
tion of these sources has been triggered by the me- 
chanical input from the first generation, 4 Myr old 
main cluster. 

6. We find a positive correlation between the strength 
of the 17 /im PAH feature, the enhancement of the 
[Sin] /[Neil] emission, and the strength of the X- 
ray field towards our sources. We propose that 



X-rays are responsible for both the excitation of 
[Sin] and the enhancement of the 17 fj,m feature 
via selective destruction of small PAH molecules. 

Our detection of molecular hydrogen in the region 
indicates gas excitation temperatures of ~ 500 K 
near selected sources, and a slightly larger abun- 
dance of H2 near field A, which correlates well 
wi th the location of a brig ht CO cloud reported 
in I Wilson fc Scovillej (fl99^ . 



The authors want to thank Mederic Boquien and the 
HERM33ES team for kindly providing the PACS pho- 
tometry of NGC 604, Dr. Inga Kamp for her comments 
on the text, and the anonymous referee for a very de- 
tailed report. This paper would not have been possible 
without the hospitality of all the astronomical commu- 
nity at the Lowell Observatory. This research has been 
possible thanks to funds provided by the Netherlands 
Research School of Astronomy (NOVA) . We made use of 
SAOImage DS9, developed by Smithsonian Astrophysi- 
cal Observatory. 



REFERENCES 



Barba, R. H., Mai'z Apellaniz, J., Perez, E., et al. 2009, Ap&SS, 
324, 309 

Beirao, P., Brandl, B. R., Devost, D., et al. 2006, ApJ, 643, LI 
Boquien, M., Calzetti, D., Combes, F., et al. 2011, AJ, 142, 111 
Brandl, B. R., et al. 2006, ApJ, 653, 1129 

Brandl, B. R., Snijders, L., den Brok, M., et al. 2009, ApJ, 699, 
1982 

Bruhweiler, F. C, Miskey, C. L., & Smith Neubig, M. 2003, AJ, 
125, 3082 

Cannon, J. M., Smith, J.-D. T., Walter, F., et al. 2006, ApJ, 647, 
293 

Cassinelli, J. P., & Swank, J. H. 1983, ApJ, 271, 681 
Chu, Y.-H., & Mac Low, M.-M. 1990, ApJ, 365, 510 
Churchwell, E., & Goss, W. M. 1999, ApJ, 514, 188 
da Cunha, E., Chariot, S., & Elbaz, D. 2008, MNRAS, 388, 1595 
Dopita, M. A., Groves, B. A., Fischera, J., et al. 2005, ApJ, 619, 
755 

Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006, ApJ, 
647, 244 

Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006, ApJS, 
167, 177 

Drissen, L., Crowther, P. A., Ubeda, L., & Martin, P. 2008, 

MNRAS, 389, 1033 
Dudik, R. P., Weingartner, J. C, Satyapal, S., et al. 2007, ApJ, 

664, 71 

Eldridge, J. J., & Relaho, M. 2011, MNRAS, 411, 235 
Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725 
Farina, C, Bosch, G. L., & Barba, R. H. 2012, AJ, 143, 43 
Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10 
Freedman, W. L., Wilson, C. D., & Madore, B. F. 1991, ApJ, 372, 
455 

Fukui, Y., & Kawamura, A. 2010, ARA&A, 48, 547 
Gonzalez Delgado, R. M., & Perez, E. 2000, MNRAS, 317, 64 
Groves, B., Dopita, M. A., Sutherland, R. S., Kewley, L. J., 

Fischera, J., Leitherer, C, Brandl, B., & van Breugel, W. 2008, 

ApJS, 176, 438 

Groves, B. A., & Allen, M. G. 2010, New Astronomy, 15, 614 
Heiner, J. S., Allen, R. J., & van der Kruit, P. C. 2009, ApJ, 700, 
545 

Helou, G., Roussel, H., Appleton, P., et al. 2004, ApJS, 154, 253 
Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 
18 

Hunter, D. A., Shaya, E. J., Holtzman, J. A., et al. 1995, ApJ, 448, 
179 

Hunter, D. A., Baum, W. A., O'Neil, E. J., Jr., & Lynds, R. 1996, 
ApJ, 456, 174 



Indebetouw, R., de Messieres, G. E., Madden, S., et al. 2009, ApJ, 
694, 84 

Kaufman, M. J., Wolfire, M. G., & Hollenbach, D. J. 2006, ApJ, 
644, 283 

Lada, C. J., & Adams, F. C. 1992, ApJ, 393, 278 
Lawton, B., Gordon, K. D., Babler, B., et al. 2010, ApJ, 716, 453 
Lcbouteiller, V., Bernard-Salas, J., Brandl, B., et al. 2008, ApJ, 
680, 398 

Levesque, E. M., Kewley, L. J., & Larson, K. L. 2010, AJ, 139, 712 
Magrini, L., Vflchez, J. M., Mampaso, A., Corradi, R. L. M., & 

Leisy, P. 2007, A&A, 470, 865 
Mai'z-Apellaniz, J. 2001, ApJ, 563, 151 

Mai'z-Apellaniz, J., Perez, E., & Mas-Hesse, J. M. 2004, AJ, 128, 
1196 

Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 
561, 766 

Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, 
ApJ, 466, 561 

Martinez-Galarza, J. R., Groves, B., Brandl, B., et al. 2011, ApJ, 

738, 176 
Melnick, J. 1980, A&A, 86, 304 

Micelotta, E. R., Jones, A. P., & Tielens, A. G. G. M. 2010, A&A, 
510, A37 

Miura, R., Okumura, S. K., Tosaki, T., et al. 2010, ApJ, 724, 1120 
O'Dowd, M. J., Schiminovich, D., Johnson, B. D., et al. 2009, ApJ, 
705, 885 

Pereira-Santaella, M., Alonso-Herrero, A., Rieke, G. H., et al. 2010, 
ApJS, 188, 447 

Poglitsch, A., Waelkens, C, Geis, N., et al. 2010, A&A, 518, L2 
Reach, W. T., Megeath, S. T., Cohen, M., et al. 2005, PASP, 117, 
978 

Relano, M., & Kennicutt, R. C. 2009, ApJ, 699, 1125 
Rubin, R. H., Simpson, J. P., Lord, S. D., et al. 1994, ApJ, 420, 
772 

Seale, J. P., Looney, L. W., Chu, Y.-H., et al. 2009, ApJ, 699, 150 
Servillat, M., Dieball, A., Webb, N. A., et al. 2008, A&A, 490, 641 
Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007, ApJ, 656, 
770 

Sofia, U. J., & Jenkins, E. B. 1998, ApJ, 499, 951 

Tenorio-Tagle, C, Munoz-Tunon, C, Perez, E., Maiz- Apellaniz, J., 

& Medina- Tanco, G. 2000, ApJ, 541, 720 
Tosaki, T., Miura, R., Sawada, T., et al. 2007, ApJ, 664, L27 
Tiillmann, R., Gaetz, T. J., Plucinsky, P. P., et al. 2008, ApJ, 685, 

919 

Van Kerckhoven, C, Hony, S., Peeters, E., et al. 2000, A&A, 357, 
1013 



Ongoing Massive Star Formation in NGC 604 



Wilson, C. D., & Scovillc, N. 1992, ApJ, 385, 512 Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 

Yang, H., Chu, Y.-H., Skillman, E. D., & Terlevich, R. 1996, AJ, 
112, 146 

Zakamska, N. L. 2010, Nature, 465, 60 



